Exercise 9.7

Problem.

Implement a Gibbs sampler to generate a bivariate normal chain \((X_t,Y_t)\) with zero means, unit standard deviations, and correlation 0.9. Plot the generated sample after discarding a suitable burn-in sample. Fit a simple linear regression model \(Y = \beta_0 + \beta_1 X\) to the sample and check the residuals of the model for normality and constant variance.

*use the Gelman-Rubin method to monitor convergence of the chain, and run the chain until it converges approximately to the target distribution according to \(\hat R < 1.2\).

Solution.

The target distribution is bivariate normal: \((X,Y)\sim N(\mu_1,\mu_2,\sigma^2_1,\sigma_2^2,\rho)\).

And the conditional distributions \(f(x|y)\) and \(f(y|x)\) are: \[(x|y)\sim N(\mu_1+\rho\sigma_1/\sigma_2(y-\mu_2),(1-\rho^2)\sigma_1^2),\] \[(y|x)\sim N(\mu_2+\rho\sigma_2/\sigma_1(x-\mu_1),(1-\rho^2)\sigma_2^2).\]

• For \(t=1,\ldots,T\)
  • 1. Sets \((x,y)=Z(t-1)\);
  • 2. Generates \(X^*(t)\) from \(f(\cdot|y)\);
  • 3. Updates \(x=X^*(t)\);
  • 4. Generates \(X^*(t)\) from \(f(\cdot|x)\);
  • 5. Sets \(X(t)=(X^*(t),Y^*(t))\).
• End for

We set length of chain: \(N=8000\),Burn-in length = 1000, \(\sigma_1^2=1\), \(\sigma^2_2=1\), \(\rho=0.9\).

First, we generate 5 bivariate normal chain using four different initial values from the codes below.

set.seed(22015)

###### generate the chain #####
bivariate_normal_chain = function(mu1,mu2,sigma1,sigma2,rho,N){
  Z = matrix(0, N, 2) #the chain, a bivariate sample
  Z[1, ] = c(mu1, mu2) #initialize
  s1 = sqrt(1-rho^2)*sigma1
  s2 = sqrt(1-rho^2)*sigma2
  for (i in 2:N) {
    y = Z[i-1, 2]
    m1 = 0 + rho * (y - 0) * sigma1/sigma2
    Z[i, 1] = rnorm(1, m1, s1)
    x = Z[i, 1]
    m2 = 0 + rho * (x - 0) * sigma2/sigma1
    Z[i, 2] = rnorm(1, m2, s2)
  }
  return(Z)
}
#initialize constants and parameters
N = 8000 #length of chain
burn = 1000 #burn-in length
rho = 0.9 #correlation
k = 5 # number of initial values
mu1 = c(-15,-10,0,5,10)
mu2 = c(-15,-10,0,5,10)
sigma1 = 1
sigma2 = 1
Z = array(0,dim = c(k,N,2))
for (i in 1:k){
  Z[i,,] = bivariate_normal_chain(mu1[i],mu2[i],sigma1,sigma2,rho,N)
}

Then I will use the Gelman-Rubin method to monitor convergence of the chains in which overdispersed initial values are chosen.

According to the General Methods for Monitoring Convergence of Iterative Simulations by Stephen P. BROOKSand Andrew GELMAN, there are a number of multivariate approaches available to us. When we are interested in the estimation of a vector parameter \(\psi\) based upon observations \(\psi_{jt}^{(i)}\) denoting the \(i\) th element of the parameter vector in chain \(j\) at time \(t\), the direct analogue of the univariate approach in higher dimensions is to estimate the posterior variance-covariance matrix by \[ \widehat{V}=\frac{n-1}{n} W+\left(1+\frac{1}{m}\right) B / n \] where \[ W=\frac{1}{m(n-1)} \sum_{j=1}^m \sum_{t=1}^n\left(\psi_{j t}-\bar{\psi}_{j .}\right)\left(\psi_{j t}-\bar{\psi}_{j .}\right)^{\prime} \]

and \[ B / n=\frac{1}{m-1} \sum_{j=1}^m\left(\bar{\psi}_{j .}-\bar{\psi}_{. .}\right)\left(\bar{\psi}_{j .}-\bar{\psi}_{. .}\right)^{\prime} \]

We would like to summarize the distance between \(\widehat{V}\) and \(W\) with a scalar measure which should approach 1 (from above) as convergence is achieved, given suitably overdispersed starting points. Many such measures can be constructed; we will work with the maximum root statistic, which is the maximum SRF of any linear projection of \(\psi\), given by \[ \widehat{R}^p=\max _a \frac{a^{\prime} \widehat{V} a}{a^{\prime} W a}, \] which we call the multivariate PSRF (or MPSRF) .

And \[ \widehat{R}^p=\frac{n-1}{n}+\left(\frac{m+1}{m}\right) \lambda_1, \] where \(\lambda_1\) is the largest eigenvalue of the symmetric, positive definite matrix \(W^{-1} B / n\).

In the article we mentioned before, Lemma 3 can help us “calculate” the multivariate PSRF.

Lemma 3. For \(k=1, \ldots, p\), let \(\widehat{R}(k)\) be the univariate PSRF , applied to the simulations of the kth element of \(\theta\). If \(\widehat{R}^{\max }\) denotes the maximum of the \(\widehat{R}(k)\) values (over \(k=1, \ldots, p)\) then, \[ \widehat{R}^{\max } \leq \widehat{R}^p \] where \(\widehat{R}^p\) is the MPSRF defined before, applied to the vector of parameters, \(\theta\).

#compute diagnostic statistics
psi = apply(Z, c(1,3), cumsum)
for (i in 1:dim(psi)[2]){
  for (j in 1:dim(psi)[3]){
    psi[,i,j] = psi[,i,j] / (1:dim(psi)[1])
  }
}

Gelman.Rubin_muti_lemma2 = function(psi){
      n = dim(psi)[1] #n
      k = dim(psi)[2] #chain numbers
      p = dim(psi)[3] #parameter numbers
      psi.w = array(0,c(p,p,k))
      for (i in 1:k){
        psi.w[,,i] = cov(psi[,i,])
      }
      
      W = apply(psi.w, c(1,2), "mean")
      B.1 = cov(apply(psi,c(2,3),"mean" ))
      #v.hat = W*(n-1)/n +B/n     
      
      # solve max singular value
       lamda = eigen(solve(W)%*%B.1)$val
       r.hat = (n-1)/n + (k+1)/k*max(lamda)
      
      return(r.hat)
}

Gelman.Rubin_muti_lemma3 = function(psi){
      n = dim(psi)[1] #n
      k = dim(psi)[2] #chain numbers
      P = dim(psi)[3] #parameter numbers
      R = rep(0,P)
      for (p in 1:P){
      data = t(psi[,,p])
      R[p] = Gelman.Rubin(data)
      }
      return(max(R)**{1/P})
}      

#plot the sequence of R-hat statistics
rhat = rep(0, N)
for (j in (burn+1):N)
    rhat[j] = Gelman.Rubin_muti_lemma3(psi[1:j,,])
plot(rhat[(burn+1):N], type="l", xlab="", ylab="R")
abline(h=1.2, lty=2)

We can see that the \(\hat R\) has been lower than 1.2 after burnin period. The scatterplot of the generated chain (after discarding the burn-in sample) is shown below.

par(mfrow=c(3,2))
plot(Z[1,(burn+1):N,], main="", cex=.5, xlab=bquote(Z[1,(burn+1):N,1]),ylab=bquote(Z[1,(burn+1):N,2]))
for (i in 2:k){
  plot(Z[i,(burn+1):N,], main="", cex=.5, xlab=bquote(Z[i,(burn+1):N,1]),ylab=bquote(Z[i,(burn+1):N,2]))
}

We fit a simple linear regression model \(Y = \beta_0 + \beta_1 X\) to the sample and check the residuals of the model for normality and constant variance.

x = c()
y = c()
for (i in 1:k){
  x = c(x,Z[i,(burn+1):N,1])
  y = c(y,Z[i,(burn+1):N,2])
}
linear_function = lm(y~x)

summary(linear_function)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.85810 -0.29288  0.00156  0.29592  1.66629 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.003077   0.002339  -1.316    0.188    
## x            0.901385   0.002328 387.135   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4373 on 34998 degrees of freedom
## Multiple R-squared:  0.8107, Adjusted R-squared:  0.8107 
## F-statistic: 1.499e+05 on 1 and 34998 DF,  p-value: < 2.2e-16

The plot of residuals vs fits suggests that the error variance is constant with respect to the response variable. The QQ plot of residuals is consistent with the normal error assumption of the linear model.

plot(linear_function)

Finally, we clean the memory of the variables.

rm(list = ls())

Back to the Question